This notebook shows how to quickly generate a slideshow from a set of cells. This can be a great way to share your ideas with others as they are being formed.
The IJulia notebook has a drop-down menu called "Cell Toolbar". You can use this toolbar to assign some slideshow function to each cell. This cell is set to "Skip" -- it will not show up when the notebook is converted to a slideshow. There are several other settings:
You can convert a notebook into slides by running the following code cell. This will generate an html file that may be rendered as a slideshow by reveal.js.
In [ ]:
# Generate slideshow from this notebook:
run(`ipython nbconvert IV.\ Exploring\ a\ new\ idea.v3.ipynb --to slides`)
Let's derive the expression for a current flow constraint.
$$\begin{align} \lvert I_{ik} \rvert &\leq I_{ik}^{max} \\ \lvert I_{ik} \rvert^2 &\leq (I_{ik}^{max})^2 \\ [Re( I_{ik} )]^2 + [Im( I_{ik} )]^2 &\leq (I_{ik}^{max})^2 \\ \end{align}$$Next we expand the real and imaginary parts of line current $I_{ik}$:
$$\begin{align} Re( I_{ik} ) &= G_{ik}[Re(E_i) - Re(E_k)] - B_{ik}[Im(E_i) - Im(E_k)] \\ Im( I_{ik} ) &= B_{ik}[Re(E_i) - Re(E_k)] + G_{ik}[Im(E_i) - Im(E_k)] \\ \end{align}$$Here we make our first assumption, zero line resistance:
$$\begin{align} Z_{ik} &= r_{ik} + jx_{ik} \\ Y_{ik} = G_{ik} + jB_{ik} &= z_{ik}^{-1} \\ r_{ik} = 0 \implies Y_{ik} &= 0 - j\frac{1}{x_{ik}} \end{align}$$Since $G_{ik}=0$ by assumption, we can simplify the expressions for real and imaginary currents:
$$\begin{align} Re( I_{ik} ) &\approxeq - B_{ik}[Im(E_i) - Im(E_k)] \\ Im( I_{ik} ) &\approxeq B_{ik}[Re(E_i) - Re(E_k)] \\ \end{align}$$Square both expressions and add them together:
$$\begin{align} Re( I_{ik} )^2 &= B_{ik}^2[Im(E_i) - Im(E_k)]^2 \\ Im( I_{ik} )^2 &= B_{ik}^2[Re(E_i) - Re(E_k)]^2 \\ \lvert I_{ik} \rvert^2 = Re( I_{ik} )^2 + Im( I_{ik} )^2 &= B_{ik}^2[\left( Im(E_i) - Im(E_k) \right)^2 + \left( Re(E_i) - Re(E_k) \right)^2] \end{align}$$Simplify:
$$\begin{align} \lvert I_{ik} \rvert^2 &= B_{ik}^2[\left( Im(E_i) - Im(E_k) \right)^2 + \left( Re(E_i) - Re(E_k) \right)^2] \\ &= B_{ik} \left[ Re(E_i)^2 + Im(E_i)^2 + Re(V_k)^2 + Im(V_k)^2 -2\left( Re(V_i)Re(V_k) + Im(V_i)Im(V_k) \right) \right] \\ &= B_{ik} \left[ V_i^2 + V_k^2 -2V_iV_k\left( \sin\theta_i\sin\theta_k + \cos\theta_i\cos\theta_k \right) \right] \\ \lvert I_{ik} \rvert^2 &= B_{ik} \left( V_i^2 + V_k^2 - 2V_iV_k\cos\theta_{ik} \right),\quad \theta_{ik} = \theta_i - \theta_k \end{align}$$Now we have an expression for the squared current magnitude on a line assuming the line's resistance is zero. How good is this assumption? Let's make a function to plot both the full and approximate current expressions on the same chart.
Since line current is a function of two voltage magnitudes, two voltage angles, and the two components of admittance ($G_{ik}$ and $B_{ik}$), there are many ways we could plot it. I have written functions to plot the following:
Let's do (1) first, and initially assume resistance is 0 (and therefore $G_{ik}=0$).
In [2]:
(Vi,Vk,Gik,Bik) = (1.1,0.9,0,-4)
plotCurrentθ(Vi,Vk,Gik,Bik);
A complete overlap. Notice that both functions are convex.
Now let's introduce non-zero resistance. A typical line has $\lvert B_{ik}/G_{ik} \rvert \approx 4$, so let's plot that case:
In [3]:
(Vi,Vk,Gik,Bik) = (1.1,0.9,1,-4)
plotCurrentθ(Vi,Vk,Gik,Bik);
The approximate function still does a good job for small angle differences; not so much for larger ones.
Now let's make some plots of type (2). Again, let's begin with resistance set to zero:
In [4]:
(θi,θk,Gik,Bik) = (0,0.2,0,-4)
plotCurrentV(θi,θk,Gik,Bik);
A perfect match. Now let's introduce some resistance:
In [5]:
(θi,θk,Gik,Bik) = (0,0.2,1,-4)
plotCurrentV(θi,θk,Gik,Bik);
Non-zero resistance adds significantly to the line current magnitude. Since the approximate line current expression neglects resistance, it tells us there is less current on the line than there actually is.
So why use an approximation that still falls short of the exact current expression? Well for one thing, the exact expression requires us to accommodate resistance throughout our analysis. It would make no sense to assume resistance is zero in the power balance equations while allowing non-zero resistance in current constraints. More importantly, though, our approximate current flow expression captures the effects of voltage magnitude differences. In contrast, the DC approximation assumes all voltages are 1 pu and therefore have no effect on line currents (since there are no angle differences to produce reactive power flow). Let's see what that looks like:
In [6]:
(θi,θk,Gik,Bik) = (0,0.2,1,-4)
plotCurrentDC(θi,θk,Gik,Bik);
The DC approximation consistently understimates line current magnitude. When the voltage magnitude difference between two buses is significant, this underestimation is severe.
What does this mean for the linear instanton problem? It means that lines may become saturated more quickly than linear instanton analysis indicates. Any line with a significant voltage magnitude difference across it is more vulnerable than linear instanton analysis tells us.
Our expression for magnitude has a cosine in it, but we can easily solve for the angle difference:
$$\begin{align} \lvert I_{ik}^2 \rvert &\approxeq \frac{V_i^2 + V_k^2 - 2V_iV_k\cos \theta_{ik}}{x_{ik}^2} \\ \implies \theta_i - \theta_k &= \cos^{-1}\left( \frac{V_i^2 + V_k^2 - (\lvert I_{ik}^{lim} \rvert x_{ik})^2}{2V_iV_k} \right) \end{align}$$Compare this constraint to the one used in the DC instanton problem:
@addConstraint(m, Congestion, θ[i] - θ[k] == sense*x[idx]*Plim[idx])
@addConstraint(m, Congestion, θ[i] - θ[k] == sense*acos((Vi^2 + Vk^2 - (Plim[idx]*x[idx])^2)/(2Vi*Vk)))
More complicated expression, but still just setting angle difference equal to constant.
Unfortunately, our expression for current magnitude has a cosine in it. This will not play well with the solver, so we need to substitute a piecewise-linear approximation for $\cos\theta_{ik}$. This is where Appendix A of Coffrin's unpublished paper comes in handy. Suppose we want to replace cosine with $s$ linear constraints. So long as the angle difference lies in the range $(-\pi/2,\pi/2)$, we can use the following procedure to obtain the following set of tangent lines:
This gives us a set of linear inequality constraints $x_{\hat{cos}}$ of the form $$\begin{align} y &\leq \sin(a)(\theta_{ik} - a) + \cos(a)~, \end{align}$$ with $a$ spaced evenly between $l$ and $h$.
In [ ]:
function solver_currentFlow_droop(psDL)
# This function uses MOSEK to perform instanton analysis for the case with conventional generator droop resopnse.
# unpack psDL (boilerplate):
(Sb,f,t,r,x,b,Y,bustype,
Gp,Gq,Dp,Dq,Rp,Rq,
Pmax,Pmin,Qmax,Qmin,Plim,
Vg,Vceiling,Vfloor,
busIdx,N,Nr,Ng,k) = unpack_psDL(psDL)
score = Float64[]
rho = Array(Vector{Float64},0) # array of vectors with Float64 values
Gpost = Array(Vector{Float64},0) # array of vectors with Float64 values
theta = Array(Vector{Float64},0) # array of vectors with Float64 values
alpha = Float64[]
result = Symbol[]
# Enforce each constraint in the f --> t direction
for idx = 1:length(f)
(m, ρ, θ, α, PowerBalance, AngleRef, Congestion, Participation) = droopCurrentModel(idx, 1, Rp, Gp, Dp, f, t, x, Y, bustype, Plim,k,Vg)
status = solve(m);
push!(score, getObjectiveValue(m))
push!(result,status)
push!(theta, getValue(θ)[:])
push!(alpha, getValue(α))
push!(rho,getValue(ρ)[:])
# Compute conventional generation post-instanton:
push!(Gpost, Gp + k.*getValue(α))
end
# Enforce each constraint in the t --> f direction
for idx = 1:length(f)
(m, ρ, θ, α, PowerBalance, AngleRef, Congestion, Participation) = droopCurrentModel(idx, -1, Rp, Gp, Dp, f, t, x, Y, bustype, Plim,k,Vg)
status = solve(m);
push!(score, getObjectiveValue(m))
push!(result, status)
push!(theta, getValue(θ)[:])
push!(alpha, getValue(α))
push!(rho, getValue(ρ)[:])
# Compute conventional generation post-instanton:
push!(Gpost, Gp + k.*getValue(α))
end
# Generate strings to tell user which constraint was violated:
constrIdx = String[]
for i = 1:2length(f)
if i <= length(f)
push!(constrIdx, "activeFlow: $(f[i]) -> $(t[i]) @ $(Plim[i])")
else
i2 = i - length(f)
push!(constrIdx, "activeFlow: $(t[i2]) -> $(f[i2]) @ $(Plim[i2])")
end
end
rankedList = sortrows([zero2NaN(score) [rho theta alpha constrIdx Gpost]]) # sort candidates; place NaNs at bottom
# Pack results into instance of "instantonResults":
return instantonResults(rankedList[:,1],rankedList[:,2],rankedList[:,3],rankedList[:,4],rankedList[:,5],rankedList[:,6])
end
function droopCurrentModel(idx, sense, Rp, Gp, Dp, f, t, x, Y, bustype, Plim,k, Vg)
# 7-29-14
# DROOP RESPONSE
# Create model saturating line 'idx' in direction 'sense' (±1)
# This function uses JuMP and Mosek
# Line current constraints
m = Model(solver = MosekSolver()) # Use MOSEK to solve model
N = length(Rp)
@defVar(m, ρ[1:N] >= 0) # Add decision variables to model (renewable gen)
@defVar(m, θ[1:N]) # Add bus angles
@defVar(m, α) # mismatch parameter
setObjective(m, :Min, 0.5*sum([(ρ[i] - Rp[i])^2 for i in find(Rp)]))
# add power balance constraints:
@defConstrRef PowerBalance[1:N]
for i in setdiff(1:N,find(bustype.==3))
PowerBalance[i] = @addConstraint(m, sum([Y[i,j]*θ[j] for j in 1:N]) == Gp[i] + k[i]*α + ρ[i] - Dp[i])
end
# wind output must be 0 at non-wind buses:
@addConstraint(m, NonWind, sum([ρ[i] for i in setdiff(collect(1:N),find(Rp))]) == 0)
# angle reference bus must have angle 0:
@addConstraint(m, AngleRef, θ[find(bustype.==3)[1]] == 0)
# 'i' is f[idx]
# 'k' is t[idx]
#i = f[idx]
#k = t[idx]
Vi = Vg[f[idx]]
Vk = Vg[t[idx]]
# saturate a line (current constraint):
@addConstraint(m, Congestion, θ[f[idx]] - θ[t[idx]] == sense*acos((Vi^2 + Vk^2 - (Plim[idx]*x[idx])^2)/(2Vi*Vk)))
#@addConstraint(m, Congestion, θ[f[idx]] - θ[t[idx]] == sense*x[idx]*Plim[idx]) # saturate a line
# @addConstraint(m, Congestion, (Vi^2 + Vk^2 - 2Vi*Vk*sense*cos(θ[i] - θ[k]))/(x[idx]^2) == Plim[idx]^2)
@addConstraint(m, Participation, α == (sum(Dp) - sum([ρ[i] for i in find(Rp)])) - sum(Gp))
return m, ρ, θ, α, PowerBalance, AngleRef, Congestion, Participation
end
Note: this part doesn't appear in the slideshow!
In [1]:
using PyPlot
PyPlot.svg(true)
function plotCurrentθ(Vi,Vk,Gik,Bik)
θi = linspace(-pi/22,pi/22,200)
θk = linspace(pi/22,-pi/22,200)
Yik = Gik - Bik*1im
Ei = [Vi*e^(θi[i]*im) for i in 1:length(θi)]
Ek = [Vk*e^(θk[i]*im) for i in 1:length(θk)]
Iik = [(Ei[i] - Ek[i])*Yik for i in 1:length(θi)]
Icalc = [(Vi^2 + Vk^2 - 2Vi*Vk*cos(θi[i] - θk[i]))*(imag(Yik)^2) for i in 1:length(θi)]
plot(θi-θk,[abs(Iik[i])^2 for i in 1:length(Iik)],linewidth=2)
hold = true
plot(θi - θk, Icalc,linewidth=2,"--")
legend(["Actual |Iik|^2","Approx. |Iik|^2"],loc=4,fontsize=10)
xlabel("Angle difference (rad)")
ylabel("Current (pu)")
xlim(-0.3,0.3)
ylim(0,3)
end
function plotCurrentV(θi,θk,Gik,Bik)
Vi = linspace(0.9,1.1,200)
Vk = linspace(1.1,0.9,200)
Yik = Gik - Bik*1im
Ei = [Vi[i]*e^(θi*im) for i in 1:length(Vi)]
Ek = [Vk[i]*e^(θk*im) for i in 1:length(Vk)]
Iik = [(Ei[i] - Ek[i])*Yik for i in 1:length(Vi)]
Icalc = [(Vi[i]^2 + Vk[i]^2 - 2Vi[i]*Vk[i]*cos(θi - θk))*(imag(Yik)^2) for i in 1:length(Vi)]
plot(Vi - Vk,[abs(Iik[i])^2 for i in 1:length(Iik)],linewidth=2)
hold = true
plot(Vi - Vk, Icalc,linewidth=2,"--")
legend(["Actual |Iik|^2","Approx. |Iik|^2"],loc=4,fontsize=10)
xlabel("Magnitude difference (pu)")
ylabel("Current (pu)")
xlim(-0.3,0.3)
ylim(0,3)
end
function plotCurrentDC(θi,θk,Gik,Bik)
Vi = linspace(0.9,1.1,200)
Vk = linspace(1.1,0.9,200)
Yik = Gik - Bik*1im
Ei = [Vi[i]*e^(θi*im) for i in 1:length(Vi)]
Ek = [Vk[i]*e^(θk*im) for i in 1:length(Vk)]
Iik = [(Ei[i] - Ek[i])*Yik for i in 1:length(Vi)]
Icalc = [((θi - θk)*imag(Yik))^2 for i in 1:length(Vi)]
plot(Vi - Vk,[abs(Iik[i])^2 for i in 1:length(Iik)],linewidth=2)
hold = true
plot(Vi - Vk, Icalc,linewidth=2,"--")
legend(["Actual |Iik|^2","DC |Iik|^2"],loc=4,fontsize=10)
xlabel("Magnitude difference (pu)")
ylabel("Current (pu)")
xlim(-0.3,0.3)
ylim(0,3)
end
Out[1]: